/* 	REALLOCATE ALL SOLAR CAPACITY
	SPLIT TENNET in N and S region
	CONSIDER TRANSMISSION CONSTRAINT 
*/


clear all
cd "${PATH}/data"


cap log close 
set logtype text
log using "../results/log_PVoutput_analysis3_splitTenneTdeltaK.txt", replace

timer clear
timer on 1
set more off
use "data_15min.dta", clear


// MERGE WITH FILE WITH LAMBDAS (AVOIDED OPERATING COSTS, dAS/dR, AVOIDED EMISSIONS) 
preserve
use "lambda_15min.dta", clear
keep TSO time_utc marginal_cost marginal_emissions kload solar step_* cost_* emission_* dAS_dR
replace time_utc = time_utc - msofhours(1)
tempfile temp1
save `temp1'
restore

merge m:1 TSO time_utc using `temp1'
keep if _merge ==3 
drop _merge 


// COMPUTE BASELINE OUTPUT
gen year = yofd(dofc(time_utc))


*** UPDATE SOLAR PRODUCTION TENNET_South AND TENNET_North
// Duplicate TenneT observations and replace "solar production data from Bavaria"
expand 2 if TSO =="TenneT", gen(dup_)
replace TSO = "TenneT_North" if TSO =="TenneT" & dup_ == 0 
replace TSO = "TenneT_South" if TSO =="TenneT" & dup_ == 1

** Data on Bavaria solar production - see file "Supply_Demand_splitTenneT.do"
preserve 
use "Tennet_solar_South.dta", clear
ren solar solar_BY
tempfile temp2
save `temp2'
restore

* Merge to all TSOs
merge m:1 time_utc using `temp2', keepusing(solar_BY)
drop if _merge == 2 // 31 december 2016
drop _merge


// UPDATE production for SOUTH / NORTH TENNET
replace solar_gen = solar_BY if TSO =="TenneT_South"
replace solar_gen = solar_gen - solar_BY if TSO =="TenneT_North"
replace solar_gen = 0 if solar_gen < 0 // very few instances when Bavaria reported more solar than TSO total
// Merging separately rather than using all from TenneT_split.dta" keeps the 15-min time intervals for solar

*******************************
// MERGE WITH THE LAMBDAS OBTAINED IN SUPPLY_DEMAND_SPLITTENNET
** Note: there are hours where supply does not intersect with demand!! ** 
gen date_hour = hh(time_utc)
gen date_day = dofc(time_utc)

merge m:1 TSO date_day date_hour using "TenneT_split.dta", keepusing(*_split lambda*)
drop if _merge == 2
drop _merge

* UPDATE "SOUTH" AND "NORTH" VALUES * 
// Use updated values for TenneT: recalculated marginal cost, marginal emissions, load
replace marginal_cost = lambda if TSO == "TenneT_North" | TSO =="TenneT_South"
replace marginal_emissions = lambda_emissions if TSO == "TenneT_North" | TSO =="TenneT_South"
replace load = load_split if TSO == "TenneT_North" | TSO =="TenneT_South"

// Update also "steps" and MC for TenneT_N and TenneT_S
forvalues s = 0/6 {
replace step_`s' = lambda_step_`s' if TSO == "TenneT_North" | TSO =="TenneT_South"
}

// Keep Marginal cost and emissions of steps for TenneT
ren lambda_mc_step_* mc_step_*
ren lambda_me_step_* me_step_*

// Need to use updated file for KLOAD CLUSTER  that contains both N and S!! 
** Update in file AS_regressions ** 
merge m:1 TSO date_day date_hour using "data15min_cluster_TSO_TenneTSplit.dta", keepusing(kload)
keep if _merge ==3 
drop _merge

encode TSO, gen(TSO_num)
xtset TSO_num time_utc, delta(15 minutes)

*********
// Need to update dAS to reflect the two TenneT areas also in the baseline!!
*********
*
drop dAS_dR

gen dAS_dR = 5.174 +  (-0.00223) * 2 * solar_gen + ///
	( -9.83e-09) * 3 * solar_gen^2  + (0.000380) * load + ///
	 ( -6.76e-08) * load^2 + ( 0.000000189) * 2 * solar_gen * load if kload == 1 & TSO == "50Hertz"
	
replace dAS_dR =  6.455 +  (-0.000391) * 2 * solar_gen + ///
	(-3.10e-08) * 3 * solar_gen^2  + (0.0000577) * load + ///
	 (-6.62e-08) * load^2 + (5.46e-08) * 2 * solar_gen * load if kload == 3 & TSO == "50Hertz"	
	
replace dAS_dR =  18.44 + (-0.00560) * 2 * solar_gen + ///
	(0.000000512) * 3 * solar_gen^2  + (-0.00260) * load + ///
	  0.000000159 * load^2 + (0.000000185) * 2 * solar_gen * load if kload == 2 & TSO == "50Hertz"			  

replace dAS_dR = 79.60 +  0.00101 * 2 * solar_gen + ///
	(0.000000101) * 3 * solar_gen^2  + (-0.00558) * load + ///
	 0.000000103 * load^2 + (-9.49e-08) * 2 * solar_gen * load if kload == 1 & TSO == "Amprion"
	
replace dAS_dR = -20.40 + (-0.000250) * 2 * solar_gen + ///
	(0.000000442) * 3 * solar_gen^2  + (0.00184) * load + ///
	 -2.90e-08 * load^2 + (-0.000000117) * 2 * solar_gen * load if kload == 2 & TSO == "Amprion"	
	
replace dAS_dR =  77.84 +  0.00388 * 2 * solar_gen + ///
	(-6.73e-09) * 3 * solar_gen^2  + (-0.00705) * load + ///
	  0.000000156 * load^2 + (-0.000000155) * 2 * solar_gen * load if kload == 3 & TSO == "Amprion"			  

replace dAS_dR =  467.9 + (-0.00219) * 2 * solar_gen + ///
	(2.35e-08) * 3 * solar_gen^2 + (-0.103) *load + ///
	(0.00000570) * load^2 + (0.000000246) *2 * solar_gen * load if kload == 1 & TSO == "TenneT_South"  

replace dAS_dR = -0.834 + (0.00221) * 2 * solar_gen + ///
	(2.13e-08) * 3 * solar_gen^2 + (-0.00168) *load + ///
	(0.000000290) * load^2 + (-0.000000367) *2 * solar_gen * load if kload == 2 & TSO == "TenneT_South"

replace dAS_dR = 66.74 + (0.00325) * 2 * solar_gen + ///
	(-0.000000162) * 3 * solar_gen^2 + (-0.0192) *load + ///
	(0.00000131) * load^2 + (-0.000000245) *2 * solar_gen * load if kload == 3 & TSO ==	"TenneT_South"

replace dAS_dR = 1052.4 + (-0.000161) * 2 * solar_gen + ///
	(0.00000105) * 3 * solar_gen^2 + ( -0.164) *load + ///
	(0.00000634) * load^2 + (-0.000000247) *2 * solar_gen * load if kload == 1 & TSO ==	 "TenneT_North" 
	  
replace dAS_dR =  9.278  + (-0.000574) * 2 * solar_gen + ///
	(-6.56e-08) * 3 * solar_gen^2 + (-0.00237) *load + ///
	(0.000000132) * load^2 + (0.000000174) *2 * solar_gen * load if kload == 2 & TSO ==  "TenneT_North" 

replace dAS_dR =  78.21 + (0.00233) * 2 * solar_gen + ///
	(-3.10e-08) * 3 * solar_gen^2 + (-0.0146) *load + ///
	(0.000000666) * load^2 + (-0.000000171) *2 * solar_gen * load if kload == 3 & TSO ==  "TenneT_North" 

replace dAS_dR = 130.8 +  0.00561 * 2 * solar_gen + ///
	(-1.23e-08) * 3 * solar_gen^2  + (-0.0329) * load + ///
	 0.00000202 * load^2 + (-0.000000616) * 2 * solar_gen * load if kload == 1 & TSO == "TransnetBW"
	
replace dAS_dR = 40.25 + 0.0160 * 2 * solar_gen + ///
	(0.00000132) * 3 * solar_gen^2  + (-0.0217) * load + ///
	 0.00000274 * load^2 + (-0.00000375) * 2 * solar_gen * load if kload == 2 & TSO == "TransnetBW"	
	
replace dAS_dR = 105.7 + 0.0130 * 2 * solar_gen + ///
	(-0.000000160) * 3 * solar_gen^2  + (-0.0314) * load + ///
	  0.00000223 * load^2 + (-0.00000150) * 2 * solar_gen * load if kload == 3 & TSO == "TransnetBW"			  
*



*******************************
// DEFINE LOCALS // 

* CHOOSE incremental step for loop, i.e. XX  MW steps reallocation at the time
local s = $step_size 

* CHOOSE GAMMA // Max share of houses to be covered with solar 
local gamma_min = $gamma_min 
local gamma_max = $gamma_max 
local gamma_delta = $gamma_delta


* CHOOSE RANGE OF DELTAS
local deltaK_start 	= $deltaK_min
local deltaK_end 	= $deltaK_max
local deltaK_inc 	= $deltaK_inc 

local number_gammas = round((`gamma_max' - `gamma_min') / `gamma_delta') + 1
local number_deltaKs = round((`deltaK_end' - `deltaK_start') / `deltaK_inc') + 1  

* CHOOSE SCC VALUE 
local SCC  $SCC 
*********************************


gen marginal_emissions_EURton = `SCC' * marginal_emissions

// Marginal costs and marginal emissions (mc, me) for each of the steps: REPLACE for other TSOs
replace mc_step_0 = 0 if regexm(TSO, "TenneT") != 1 // Renewables (zero MC)
replace mc_step_1 = 5.1 if regexm(TSO, "TenneT") != 1 // Nuclear 
replace mc_step_2 = 6 if regexm(TSO, "TenneT") != 1 // Waste and Biomass 
replace mc_step_3 = 11.3  if regexm(TSO, "TenneT") != 1 // Lignite
replace mc_step_4 = cost_MWH_coal if regexm(TSO, "TenneT") != 1
replace mc_step_5 = cost_MWH_ng if regexm(TSO, "TenneT") != 1
replace mc_step_6 = cost_MWH_crude_oil if regexm(TSO, "TenneT") != 1

replace me_step_0 = 0 if regexm(TSO, "TenneT") != 1 // Renewables (zero MC + zero emission)
replace me_step_1 = 0 if regexm(TSO, "TenneT") != 1 // Nuclear (no emissions)
replace me_step_2 = 0 if regexm(TSO, "TenneT") != 1 // Waste and Biomass (no emissions)
replace me_step_3 = 1.148*`SCC' if regexm(TSO, "TenneT") != 1 // Lignite if regexm(TSO, "TenneT") != 1
replace me_step_4 = emission_MWH_coal*`SCC' if regexm(TSO, "TenneT") != 1 // hard coal
replace me_step_5 = emission_MWH_ng*`SCC' if regexm(TSO, "TenneT") != 1 // gas
replace me_step_6 = emission_MWH_crude_oil*`SCC' if regexm(TSO, "TenneT") != 1 // oil

drop cost_* emission_* solar

// DEFINE GAMMA 
** GAMMA AS % OF TOTAL CAPACITY ** 

// Total installed capacities
local K_50Hertz 50082
local K_Amprion 66049
local K_TransnetBW 22096
local K_North  	44135 
local K_South  	20442

gen 	max_cap = `K_50Hertz' if TSO == "50Hertz" 
replace max_cap = `K_Amprion' if TSO == "Amprion" 
replace max_cap = `K_TransnetBW' if TSO == "TransnetBW" 
replace max_cap = `K_North' if TSO == "TenneT_North" 
replace max_cap = `K_South' if TSO == "TenneT_South" 

gen max_cap_gamma = . // max_cap * gamma (in loop)


// LIMIT sample only to positive solar output hours
keep if solar_gen >0 & !missing(solar_gen) 
****************

// Production of residential installations in MWh
rename solar_gen solar_gen_official
gen solar_gen  = .
gen solar_gen_MW = .


local solarcap_all_50Hertz 		9842.2686
local solarcap_all_Amprion  	9549.0361
local solarcap_all_TenneT_North	9058.0381 
local solarcap_all_TenneT_South	6454.5093
local solarcap_all_TransnetBW	5605.6885

// S IS THE SUM OF THE solarcap_10kw_*
local S  = 40509 // 40509.541 // SUM solarcarp_all (max values 2016)


levelsof TSO, local(T)
foreach t in `T' {
replace solar_gen_MW = solar_gen_official / `solarcap_all_`t'' if TSO =="`t'" // production of res. systems [per MW]
}

// For reallocation: block s of solar capacity (MW) that is reallocated // 
gen solar_gen_s = solar_gen_MW * `s'

*********

// Calculate Benchmark value: total and by TSO
** Values from lambda_v1.do (lambda.dta)
gen MB_bench  = -dAS_dR + marginal_cost +  marginal_emissions_EURton

*Calculate total value: MB*solar output
gen MB_bench_tot = MB_bench * solar_gen_official
gen MB_bench_AS = -dAS_dR * solar_gen_official
gen MB_bench_MC = marginal_cost * solar_gen_official
gen MB_bench_ME = marginal_emissions_EURton * solar_gen_official


// AGGREGATE TOTAL BENEFITS 
preserve

	// Collapse by TSO
	collapse (sum) MB_bench_* , by(TSO)
	
	qui: su MB_bench_tot  if TSO == "50Hertz"
		local value_50Hertz_bench = r(mean)
	qui: su  MB_bench_tot   if TSO == "Amprion"	
		local value_Amprion_bench = r(mean)
	qui: su  MB_bench_tot   if TSO == "TenneT_North"	
		local value_TenneT_N_bench = r(mean)
	qui: su  MB_bench_tot   if TSO == "TenneT_South"	
		local value_TenneT_S_bench = r(mean)
	qui: su  MB_bench_tot   if TSO == "TransnetBW"
		local value_TransnetBW_bench = r(mean)
 
	qui: su MB_bench_AS  if TSO == "50Hertz"
		local MB_bench_AS_50Hertz = r(mean)
	qui: su MB_bench_AS  if TSO == "Amprion"
		local MB_bench_AS_Amprion = r(mean)
	qui: su MB_bench_AS  if TSO == "TenneT_North"
		local MB_bench_AS_TenneT_N = r(mean)
	qui: su MB_bench_AS  if TSO == "TenneT_South"
		local MB_bench_AS_TenneT_S = r(mean)
	qui: su MB_bench_AS  if TSO == "TransnetBW"
		local MB_bench_AS_TransnetBW = r(mean)

	qui: su MB_bench_MC  if TSO == "50Hertz"
		local MB_bench_MC_50Hertz = r(mean)
	qui: su MB_bench_MC  if TSO == "Amprion"
		local MB_bench_MC_Amprion = r(mean)
	qui: su MB_bench_MC  if TSO == "TenneT_North"
		local MB_bench_MC_TenneT_N = r(mean)
	qui: su MB_bench_MC  if TSO == "TenneT_South"
		local MB_bench_MC_TenneT_S = r(mean)
	qui: su MB_bench_MC  if TSO == "TransnetBW"
		local MB_bench_MC_TransnetBW = r(mean)
		
	qui: su MB_bench_ME  if TSO == "50Hertz"
		local MB_bench_ME_50Hertz = r(mean)
	qui: su MB_bench_ME  if TSO == "Amprion"
		local MB_bench_ME_Amprion = r(mean)
	qui: su MB_bench_ME  if TSO == "TenneT_North"
		local MB_bench_ME_TenneT_N = r(mean)
	qui: su MB_bench_ME  if TSO == "TenneT_South"
		local MB_bench_ME_TenneT_S = r(mean)
	qui: su MB_bench_ME  if TSO == "TransnetBW"
		local MB_bench_ME_TransnetBW = r(mean)

	
	// Collapse for total value
	collapse (sum) MB_bench_* 
	qui: su MB_bench_tot
		local value_bench_total = r(mean)
	qui: su MB_bench_AS
		local value_bench_AS = r(mean)
	qui: su MB_bench_MC
		local value_bench_MC = r(mean)
	qui: su MB_bench_ME
		local value_bench_ME = r(mean)
	
restore

drop MB_bench MB_bench_* dAS_dR 

*********
// Generate Matrices to store values
*Matrix 1: Allocation, Value. by TSO, total
mat ratios = J(`number_gammas'*`number_deltaKs', 25,.)

*********
// Reallocation //

** SAVE MAIN FILE BEFORE K loop
cap erase temp.dta

compress
save temp.dta

local m 1

// LOOP FOR TRANSMISSION CAPACITY
forvalues deltaK = `deltaK_start'(`deltaK_inc')`deltaK_end' {
di " " 
di " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% " 
di " Current level of delta K: `deltaK'"

// Start with new data for each iteration deltaK
use temp.dta, clear

local j 1

//LOOP FOR CAPACITY SHARE OF SOLAR
forvalues gamma = `gamma_min'(`gamma_delta')`gamma_max' { 	// LOOP OVER GAMMA VALUES
		di " " 
		di " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% " 
		di " Current level of gamma: `gamma' "

	*Generate variable to keep track of installed SOLAR in reallocation  
	cap drop solar_cap solar_capl1  cum_solar_gen solar_new value_50Hertz  value_Amprion value_TenneT_North value_TenneT_South value_TransnetBW
	cap drop value_AS_* value_MC_* value_ME_*
	gen solar_cap = 0 // solar capacity in loop (i)
	gen solar_capl1 = 0 // solar capacity at iteration (i-1)
	gen cum_solar_gen = 0 // keep track of cumulative solar production -> for correct "step" on supply curve
	gen solar_new = . // define difference of solar_reallocated - solar_original 
	gen lambda_S1 = . // Store position of lambda_S in loop 
	
	// Store values: Total and decomposition: ASC, MC, ME
	gen value_50Hertz = 0
	gen value_Amprion = 0
	gen value_TenneT_North = 0
	gen value_TenneT_South = 0
	gen value_TransnetBW = 0
	
	gen value_AS_50Hertz = 0
	gen value_AS_Amprion = 0
	gen value_AS_TenneT_North = 0
	gen value_AS_TenneT_South = 0
	gen value_AS_TransnetBW = 0
	gen value_MC_50Hertz = 0
	gen value_MC_Amprion = 0
	gen value_MC_TenneT_North = 0
	gen value_MC_TenneT_South = 0
	gen value_MC_TransnetBW = 0
	gen value_ME_50Hertz = 0
	gen value_ME_Amprion = 0
	gen value_ME_TenneT_North = 0
	gen value_ME_TenneT_South = 0
	gen value_ME_TransnetBW = 0
	
*************************************************	
	
	// Set max. capacity per TSO 
	replace max_cap_gamma = max_cap*`gamma'

	// LOOP OVER MW Blocks
	forvalues i = `s'(`s')`S' {  // loop over all MW of installed capacity to determine optimal allocation
		di " " 
		di " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% " 
		di "Reallocating MW: `i' of `S'"

	
		// Cumulative solar generation: include additional MW of solar installed
		replace solar_cap = solar_cap + `s'
		replace cum_solar_gen =  solar_gen_MW * solar_cap

		
		** Ancillary Costs: for last installed MW
		** RELABEL solar_gen 'cum_solar_gen' **

		gen dAS_dR = 5.174 +  (-0.00223) * 2 * cum_solar_gen + ///
			( -9.83e-09) * 3 * cum_solar_gen^2  + (0.000380) * load + ///
			 ( -6.76e-08) * load^2 + ( 0.000000189) * 2 * cum_solar_gen * load if kload == 1 & TSO == "50Hertz"
			
		replace dAS_dR =  6.455 +  (-0.000391) * 2 * cum_solar_gen + ///
			(-3.10e-08) * 3 * cum_solar_gen^2  + (0.0000577) * load + ///
			 (-6.62e-08) * load^2 + (5.46e-08) * 2 * cum_solar_gen * load if kload == 3 & TSO == "50Hertz"	
			
		replace dAS_dR =  18.44 + (-0.00560) * 2 * cum_solar_gen + ///
			(0.000000512) * 3 * cum_solar_gen^2  + (-0.00260) * load + ///
			  0.000000159 * load^2 + (0.000000185) * 2 * cum_solar_gen * load if kload == 2 & TSO == "50Hertz"			  

		replace dAS_dR = 79.60 +  0.00101 * 2 * cum_solar_gen + ///
			(0.000000101) * 3 * cum_solar_gen^2  + (-0.00558) * load + ///
			 0.000000103 * load^2 + (-9.49e-08) * 2 * cum_solar_gen * load if kload == 1 & TSO == "Amprion"
			
		replace dAS_dR = -20.40 + (-0.000250) * 2 * cum_solar_gen + ///
			(0.000000442) * 3 * cum_solar_gen^2  + (0.00184) * load + ///
			 -2.90e-08 * load^2 + (-0.000000117) * 2 * cum_solar_gen * load if kload == 2 & TSO == "Amprion"	
			
		replace dAS_dR =  77.84 +  0.00388 * 2 * cum_solar_gen + ///
			(-6.73e-09) * 3 * cum_solar_gen^2  + (-0.00705) * load + ///
			  0.000000156 * load^2 + (-0.000000155) * 2 * cum_solar_gen * load if kload == 3 & TSO == "Amprion"			  

		replace dAS_dR =  467.9 + (-0.00219) * 2 * cum_solar_gen + ///
			(2.35e-08) * 3 * cum_solar_gen^2 + (-0.103) *load + ///
			(0.00000570) * load^2 + (0.000000246) *2 * cum_solar_gen * load if kload == 1 & TSO == "TenneT_South"  

		replace dAS_dR = -0.834 + (0.00221) * 2 * cum_solar_gen + ///
			(2.13e-08) * 3 * cum_solar_gen^2 + (-0.00168) *load + ///
			(0.000000290) * load^2 + (-0.000000367) *2 * cum_solar_gen * load if kload == 2 & TSO == "TenneT_South"

		replace dAS_dR = 66.74 + (0.00325) * 2 * cum_solar_gen + ///
			(-0.000000162) * 3 * cum_solar_gen^2 + (-0.0192) *load + ///
			(0.00000131) * load^2 + (-0.000000245) *2 * cum_solar_gen * load if kload == 3 & TSO ==	"TenneT_South"

		replace dAS_dR = 1052.4 + (-0.000161) * 2 * cum_solar_gen + ///
			(0.00000105) * 3 * cum_solar_gen^2 + ( -0.164) *load + ///
			(0.00000634) * load^2 + (-0.000000247) *2 * cum_solar_gen * load if kload == 1 & TSO ==	 "TenneT_North" 
			  
		replace dAS_dR =  9.278  + (-0.000574) * 2 * cum_solar_gen + ///
			(-6.56e-08) * 3 * cum_solar_gen^2 + (-0.00237) *load + ///
			(0.000000132) * load^2 + (0.000000174) *2 * cum_solar_gen * load if kload == 2 & TSO ==  "TenneT_North" 

		replace dAS_dR =  78.21 + (0.00233) * 2 * cum_solar_gen + ///
			(-3.10e-08) * 3 * cum_solar_gen^2 + (-0.0146) *load + ///
			(0.000000666) * load^2 + (-0.000000171) *2 * cum_solar_gen * load if kload == 3 & TSO ==  "TenneT_North" 

		replace dAS_dR = 130.8 +  0.00561 * 2 * cum_solar_gen + ///
			(-1.23e-08) * 3 * cum_solar_gen^2  + (-0.0329) * load + ///
			 0.00000202 * load^2 + (-0.000000616) * 2 * cum_solar_gen * load if kload == 1 & TSO == "TransnetBW"
			
		replace dAS_dR = 40.25 + 0.0160 * 2 * cum_solar_gen + ///
			(0.00000132) * 3 * cum_solar_gen^2  + (-0.0217) * load + ///
			 0.00000274 * load^2 + (-0.00000375) * 2 * cum_solar_gen * load if kload == 2 & TSO == "TransnetBW"	
			
		replace dAS_dR = 105.7 + 0.0130 * 2 * cum_solar_gen + ///
			(-0.000000160) * 3 * cum_solar_gen^2  + (-0.0314) * load + ///
			  0.00000223 * load^2 + (-0.00000150) * 2 * cum_solar_gen * load if kload == 3 & TSO == "TransnetBW"			  

			
		// UPDATE SECTION ON INTERCONNECTION: need to do within loop to identify correct lambda_S and lambda_N  // 
		** IMPUTE VALUES FOR LAMBDA_N AS A FUNCTION OF THE CAPACITY CONSTRAINT AND LAMBDA_S
		// From file Supply_Demand_splitTenneT.do use estimates with gap = 2
		local coeff1    -0.000932 
		local coeff2    -0.00634 
			
		// Define location on supply curve -> identify MC and marginal Emissions
		replace solar_new =  cum_solar_gen - solar_gen_official // how much more solar do I generate compared to original allocation? 
		** SUBSTITUTE WITH SOLAR GEN_OFFICIAL HERE
		
		
		if (solar_new <= 0) { 
			*Apply E[MB] as "observed": assumption of "flat MB curves"
			*di "%%%%%%%%%%%%%%%%%%%%%%%%%%"
			*di "In flat part of MB"
			gen MB = -dAS_dR + marginal_cost + marginal_emissions_EURton
			
			** Assign value to production: all TSOs
			gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) // evaluate MB for block of reallocated solar
			gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0)
			gen MB_MC = cond(cum_solar_gen  < load, marginal_cost * solar_gen_s, 0)
			gen MB_ME = cond(cum_solar_gen  < load, marginal_emissions_EURton * solar_gen_s , 0)
	
			** Generate variable to keep track of solar_new_South and carryforward to all TSO 
			gen solar_new_South = solar_new if TSO == "TenneT_South"
			gsort time_utc -solar_new_South
			by time_utc: replace solar_new_South = solar_new_South[_n-1] if missing(solar_new_South) // within same time interval, carryforward to all TSOs

			
		*** SECTION ON SIMULATED LAMBDA FOR NORTH - TRANSMISSION **** 
			// Calculate lambda_N_simulated (fixed for given s, K, gamma)
			gen z_t = marginal_cost - lambda_S // use lagged value in expression for lambda_N_sim
			sort TSO_num time_utc
			gen lag_z_t = l.z_t			
			gen lambda_N_simulated = lambda_S + lag_z_t + (`coeff1' - `coeff2') * `deltaK' if TSO == "TenneT_North"	

			//Depending if excess + transfer to North, adjust 
			gen MB_North_sim = -dAS_dR + lambda_N_simulated + marginal_emissions_EURton if TSO == "TenneT_North"
			
			// Carryforward this variable, so we can use it in calculation for solar_value in South
			gsort time_utc -MB_North_sim
			by time_utc: replace MB_North_sim = MB_North_sim[_n-1] if missing(MB_North_sim)
			
			// Determine possible energy excess in South and adjust value
			gen excess_south = cum_solar_gen - load if  TSO == "TenneT_South"
			gen d_excess_south = excess_south > 0 & !missing(excess_south)
			replace excess_south = 0 if missing(excess_south) 
			
		* CONDITION 1: EXCESS ENERGY + CAPACITY CONSTRAINT NOT BINDING
			gen cond = (excess_south > 0 & excess_south <= `deltaK') //  only TenneT_S
			gsort time_utc -cond
			by time_utc: replace cond = cond[_n-1] if missing(cond) // within same time interval, carryforward to all TSOs

			*Adjust values for South (export) and North (apply same MB_sim)
			replace MB_tot = excess_south * MB_North_sim if TSO == "TenneT_South" & cond == 1	
			replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond == 1	

			
		* CONDITION 2: EXCESS ENERGY + CAPACITY CONSTRAINT BINDING
			gen cond_congest = (excess_south > 0 & excess_south > `deltaK') // constructed only for TenneT_S
			gsort time_utc -cond_congest
			by time_utc: replace cond_congest = cond_congest[_n-1] if missing(cond_congest) // within same time interval, carryforward to all TSOs
	
			*Adjust values for South (export) and North (apply same MB_sim)
			replace MB_tot = `deltaK' * MB_North_sim if TSO == "TenneT_South" & cond_congest == 1	
			replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond_congest == 1	
			
		*** END SECTION ON TRANSMISSION 
		}
	
	
		else {
			*We reallocate more solar to TSO as in "benchmark" case -> need to adjust benefits
	
			if ((solar_new > 0) & (solar_new <= step_6)) {
				* Apply MB from step6
				* di "In step 6"

				gen MB = -dAS_dR + mc_step_6 + me_step_6
				
				gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
				gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0)
				gen MB_MC = cond(cum_solar_gen  < load, mc_step_6 * solar_gen_s, 0)
				gen MB_ME = cond(cum_solar_gen  < load, me_step_6 * solar_gen_s, 0)
				* local in_step_6 `in_step_6' + 1
				
				
			*** SECTION ON SIMULATED LAMBDA FOR NORTH - TRANSMISSION **** 
				// NEW: NEED TO DEFINE WHICH lambda_S applies!!
				* While lambda_N = mc6, lambda_S can take on any value from MC0 to MC6: if-statements * 
				
				** Generate variable to keep track of solar_new_South and carryforward to all TSO 
				gen solar_new_South = solar_new if TSO == "TenneT_South"
				gsort time_utc -solar_new_South
				by time_utc: replace solar_new_South = solar_new_South[_n-1] if missing(solar_new_South) // within same time interval, carryforward to all TSOs
				
				if ((solar_new_South > 0) & (solar_new_South <= lambda_S_step_6)) {
				replace lambda_S1 = lambda_S_step_6 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6) & (solar_new_South <= (lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_5 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 ) & (solar_new_South <= (lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_4 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4  ) & (solar_new_South <= (lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_3 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 ) & (solar_new_South <= (lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_2 if TSO == "TenneT_North" 
				}				

				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 + lambda_S_step_2 ) & (solar_new_South <= (lambda_S_step_1 + lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_1 if TSO == "TenneT_North" 
				}
				
				else {
				replace lambda_S1 = lambda_S_step_0 if TSO == "TenneT_North" 
				}
				
				gen z_t = mc_step_6 - lambda_S1 // defined as lambda_N - lambda_S
				sort TSO_num time_utc
				gen lag_z_t = l.z_t			
				
				gen lambda_N_simulated = lambda_S1 + lag_z_t + (`coeff1' - `coeff2') * `deltaK' if TSO == "TenneT_North"	
				gen MB_North_sim = -dAS_dR + lambda_N_simulated + marginal_emissions_EURton if TSO == "TenneT_North"
				
				// Carryforward this variable, so we can use it in calculation for solar_value in South
				gsort time_utc -MB_North_sim
				by time_utc: replace MB_North_sim = MB_North_sim[_n-1] if missing(MB_North_sim)
				
				// Determine possible energy excess in South and adjust value
				gen excess_south = cum_solar_gen - load if  TSO == "TenneT_South"
				gen d_excess_south = excess_south > 0 & !missing(excess_south)
				replace excess_south = 0 if missing(excess_south) 
				
			* CONDITION 1: EXCESS ENERGY + CAPACITY CONSTRAINT NOT BINDING
				gen cond = (excess_south > 0 & excess_south <= `deltaK') 
				gsort time_utc -cond
				by time_utc: replace cond = cond[_n-1] if missing(cond) // within same time interval, carryforward to all TSOs

				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = excess_south * MB_North_sim if TSO == "TenneT_South" & cond == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond == 1	

				
			* CONDITION 2: EXCESS ENERGY + CAPACITY CONSTRAINT BINDING
				gen cond_congest = (excess_south > 0 & excess_south > `deltaK') 
				gsort time_utc -cond_congest
				by time_utc: replace cond_congest = cond_congest[_n-1] if missing(cond_congest) // within same time interval, carryforward to all TSOs
		
				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = `deltaK' * MB_North_sim if TSO == "TenneT_South" & cond_congest == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond_congest == 1	
				
			*** END SECTION ON TRANSMISSION 		
			}
	
	
			else if ((solar_new > step_6) & (solar_new <= (step_5 + step_6))) {
				*Apply MB from step 5
				*di "%%%%%%%%%%%%%%%%%%%%%%%%%%"
				*di "In step 5"

				gen MB = -dAS_dR + mc_step_5 + me_step_5
				gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
				gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0 )
				gen MB_MC = cond(cum_solar_gen  < load, mc_step_5 * solar_gen_s, 0)
				gen MB_ME = cond(cum_solar_gen  < load, me_step_5 * solar_gen_s, 0)
				*local in_step_5 `in_step_5' + 1
				
			*** SECTION ON SIMULATED LAMBDA FOR NORTH - TRANSMISSION **** 
				// NEW: NEED TO DEFINE WHICH lambda_S applies!!
				* While lambda_N = mc6, lambda_S can take on any value from MC0 to MC6: check in loop! * 
				
				** Generate variable to keep track of solar_new_South and carryforward to all TSO 
				gen solar_new_South = solar_new if TSO == "TenneT_South"
				gsort time_utc -solar_new_South
				by time_utc: replace solar_new_South = solar_new_South[_n-1] if missing(solar_new_South) // within same time interval, carryforward to all TSOs
				
				if ((solar_new_South > 0) & (solar_new_South <= lambda_S_step_6)) {
				replace lambda_S1 = lambda_S_step_6 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6) & (solar_new_South <= (lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_5 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 ) & (solar_new_South <= (lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_4 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4  ) & (solar_new_South <= (lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_3 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 ) & (solar_new_South <= (lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_2 if TSO == "TenneT_North" 
				}				

				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 + lambda_S_step_2 ) & (solar_new_South <= (lambda_S_step_1 + lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_1 if TSO == "TenneT_North" 
				}
				
				else {
				replace lambda_S1 = lambda_S_step_0 if TSO == "TenneT_North" 
				}
				
				gen z_t = mc_step_5 - lambda_S1 
				sort TSO_num time_utc				
				gen lag_z_t = l.z_t			
				
				gen lambda_N_simulated = lambda_S1 + lag_z_t + (`coeff1' - `coeff2') * `deltaK' if TSO == "TenneT_North"	
				gen MB_North_sim = -dAS_dR + lambda_N_simulated + marginal_emissions_EURton if TSO == "TenneT_North"
				
				// Carryforward this variable, so we can use it in calculation for solar_value in South
				gsort time_utc -MB_North_sim
				by time_utc: replace MB_North_sim = MB_North_sim[_n-1] if missing(MB_North_sim)
				
				// Determine possible energy excess in South and adjust value
				gen excess_south = cum_solar_gen - load if  TSO == "TenneT_South"
				gen d_excess_south = excess_south > 0 & !missing(excess_south)
				replace excess_south = 0 if missing(excess_south) 
				
			* CONDITION 1: EXCESS ENERGY + CAPACITY CONSTRAINT NOT BINDING
				gen cond = (excess_south > 0 & excess_south <= `deltaK') 
				gsort time_utc -cond
				by time_utc: replace cond = cond[_n-1] if missing(cond) // within same time interval, carryforward to all TSOs

				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = excess_south * MB_North_sim if TSO == "TenneT_South" & cond == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond == 1	

				
			* CONDITION 2: EXCESS ENERGY + CAPACITY CONSTRAINT BINDING
				gen cond_congest = (excess_south > 0 & excess_south > `deltaK') 
				gsort time_utc -cond_congest
				by time_utc: replace cond_congest = cond_congest[_n-1] if missing(cond_congest) // within same time interval, carryforward to all TSOs
		
				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = `deltaK' * MB_North_sim if TSO == "TenneT_South" & cond_congest == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond_congest == 1	
				
			*** END SECTION ON TRANSMISSION 	

			}
	
	
			else if  ((solar_new > step_6 + step_5) & (solar_new <= (step_4 + step_5 + step_6))) {
				* Apply MB from step 4
				*di "%%%%%%%%%%%%%%%%%%%%%%%%%%"
				*di "In step 4"
				gen MB = -dAS_dR + mc_step_4 + me_step_4
				gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
				gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0 )
				gen MB_MC = cond(cum_solar_gen  < load, mc_step_4 * solar_gen_s, 0)
				gen MB_ME = cond(cum_solar_gen  < load, me_step_4 * solar_gen_s, 0)
				*local in_step_4 `in_step_4' + 1
				
			*** SECTION ON SIMULATED LAMBDA FOR NORTH - TRANSMISSION **** 
				// NEW: NEED TO DEFINE WHICH lambda_S applies!!
				* While lambda_N = mc6, lambda_S can take on any value from MC0 to MC6: check in loop! * 
				
				** Generate variable to keep track of solar_new_South and carryforward to all TSO 
				gen solar_new_South = solar_new if TSO == "TenneT_South"
				gsort time_utc -solar_new_South
				by time_utc: replace solar_new_South = solar_new_South[_n-1] if missing(solar_new_South) // within same time interval, carryforward to all TSOs
				
				if ((solar_new_South > 0) & (solar_new_South <= lambda_S_step_6)) {
				replace lambda_S1 = lambda_S_step_6 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6) & (solar_new_South <= (lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_5 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5) & (solar_new_South <= (lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_4 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4  ) & (solar_new_South <= (lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_3 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 ) & (solar_new_South <= (lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_2 if TSO == "TenneT_North" 
				}				

				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 + lambda_S_step_2 ) & (solar_new_South <= (lambda_S_step_1 + lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_1 if TSO == "TenneT_North" 
				}
				
				else {
				replace lambda_S1 = lambda_S_step_0 if TSO == "TenneT_North" 
				}
				
				gen z_t = mc_step_4 - lambda_S1 
				sort TSO_num time_utc
				gen lag_z_t = l.z_t			
				
				gen lambda_N_simulated = lambda_S1 + lag_z_t + (`coeff1' - `coeff2') * `deltaK' if TSO == "TenneT_North"	
				gen MB_North_sim = -dAS_dR + lambda_N_simulated + marginal_emissions_EURton if TSO == "TenneT_North"
				
				// Carryforward this variable, so we can use it in calculation for solar_value in South
				gsort time_utc -MB_North_sim
				by time_utc: replace MB_North_sim = MB_North_sim[_n-1] if missing(MB_North_sim)
				
				// Determine possible energy excess in South and adjust value
				gen excess_south = cum_solar_gen - load if  TSO == "TenneT_South"
				gen d_excess_south = excess_south > 0 & !missing(excess_south)
				replace excess_south = 0 if missing(excess_south) 
				
			* CONDITION 1: EXCESS ENERGY + CAPACITY CONSTRAINT NOT BINDING
				gen cond = (excess_south > 0 & excess_south <= `deltaK') 
				gsort time_utc -cond
				by time_utc: replace cond = cond[_n-1] if missing(cond) // within same time interval, carryforward to all TSOs

				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = excess_south * MB_North_sim if TSO == "TenneT_South" & cond == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond == 1	

				
			* CONDITION 2: EXCESS ENERGY + CAPACITY CONSTRAINT BINDING
				gen cond_congest = (excess_south > 0 & excess_south > `deltaK') 
				gsort time_utc -cond_congest
				by time_utc: replace cond_congest = cond_congest[_n-1] if missing(cond_congest) // within same time interval, carryforward to all TSOs
		
				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = `deltaK' * MB_North_sim if TSO == "TenneT_South" & cond_congest == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond_congest == 1	
				
			*** END SECTION ON TRANSMISSION 	
			}
	
	
			else if  ((solar_new > step_6 + step_5 + step_4) & (solar_new <= (step_3 + step_4 + step_5 + step_6))) {
				* Apply MB from step 3
				*di "%%%%%%%%%%%%%%%%%%%%%%%%%%"
				*di "In step 3"
				gen MB = -dAS_dR + mc_step_3 + me_step_3
				gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
				gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0)
				gen MB_MC = cond(cum_solar_gen  < load, mc_step_3 * solar_gen_s, 0)
				gen MB_ME = cond(cum_solar_gen  < load, me_step_3 * solar_gen_s, 0)
				*local in_step_3 `in_step_3' + 1
				
			*** SECTION ON SIMULATED LAMBDA FOR NORTH - TRANSMISSION **** 
				// NEW: NEED TO DEFINE WHICH lambda_S applies!!
				* While lambda_N = mc6, lambda_S can take on any value from MC0 to MC6: check in loop! * 
				
				** Generate variable to keep track of solar_new_South and carryforward to all TSO 
				gen solar_new_South = solar_new if TSO == "TenneT_South"
				gsort time_utc -solar_new_South
				by time_utc: replace solar_new_South = solar_new_South[_n-1] if missing(solar_new_South) // within same time interval, carryforward to all TSOs
				
				if ((solar_new_South > 0) & (solar_new_South <= lambda_S_step_6)) {
				replace lambda_S1 = lambda_S_step_6 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6) & (solar_new_South <= (lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_5 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 ) & (solar_new_South <= (lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_4 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4  ) & (solar_new_South <= (lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_3 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 ) & (solar_new_South <= (lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_2 if TSO == "TenneT_North" 
				}				

				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 + lambda_S_step_2 ) & (solar_new_South <= (lambda_S_step_1 + lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_1 if TSO == "TenneT_North" 
				}
				
				else {
				replace lambda_S1 = lambda_S_step_0 if TSO == "TenneT_North" 
				}
				
				gen z_t = mc_step_3 - lambda_S1 
				sort TSO_num time_utc
				gen lag_z_t = l.z_t			
				
				gen lambda_N_simulated = lambda_S1 + lag_z_t + (`coeff1' - `coeff2') * `deltaK' if TSO == "TenneT_North"	
				gen MB_North_sim = -dAS_dR + lambda_N_simulated + marginal_emissions_EURton if TSO == "TenneT_North"
				
				// Carryforward this variable, so we can use it in calculation for solar_value in South
				gsort time_utc -MB_North_sim
				by time_utc: replace MB_North_sim = MB_North_sim[_n-1] if missing(MB_North_sim)
				
				// Determine possible energy excess in South and adjust value
				gen excess_south = cum_solar_gen - load if  TSO == "TenneT_South"
				gen d_excess_south = excess_south > 0 & !missing(excess_south)
				replace excess_south = 0 if missing(excess_south) 
				
			* CONDITION 1: EXCESS ENERGY + CAPACITY CONSTRAINT NOT BINDING
				gen cond = (excess_south > 0 & excess_south <= `deltaK')
				gsort time_utc -cond
				by time_utc: replace cond = cond[_n-1] if missing(cond) // within same time interval, carryforward to all TSOs

				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = excess_south * MB_North_sim if TSO == "TenneT_South" & cond == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond == 1	
				
			* CONDITION 2: EXCESS ENERGY + CAPACITY CONSTRAINT BINDING
				gen cond_congest = (excess_south > 0 & excess_south > `deltaK')
				gsort time_utc -cond_congest
				by time_utc: replace cond_congest = cond_congest[_n-1] if missing(cond_congest) // within same time interval, carryforward to all TSOs
		
				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = `deltaK' * MB_North_sim if TSO == "TenneT_South" & cond_congest == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond_congest == 1	
				
			*** END SECTION ON TRANSMISSION 		

			}
	
			else if  ((solar_new > step_6 + step_5 + step_4 + step_3) & (solar_new <= (step_2 + step_3 + step_4 + step_5 + step_6))) {
				* Apply MB from step 2
				*di "%%%%%%%%%%%%%%%%%%%%%%%%%%"
				*di "In step 2"
				gen MB = -dAS_dR + mc_step_2 + me_step_2
				gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
				gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0 )
				gen MB_MC = cond(cum_solar_gen  < load, mc_step_2 * solar_gen_s, 0)
				gen MB_ME = cond(cum_solar_gen  < load, me_step_2 * solar_gen_s, 0)
				*local in_step_2 `in_step_2' + 1
				
			*** SECTION ON SIMULATED LAMBDA FOR NORTH - TRANSMISSION **** 
				// NEW: NEED TO DEFINE WHICH lambda_S applies!!
				* While lambda_N = mc6, lambda_S can take on any value from MC0 to MC6: check in loop! * 
				
				** Generate variable to keep track of solar_new_South and carryforward to all TSO 
				gen solar_new_South = solar_new if TSO == "TenneT_South"
				gsort time_utc -solar_new_South
				by time_utc: replace solar_new_South = solar_new_South[_n-1] if missing(solar_new_South) // within same time interval, carryforward to all TSOs
				
				if ((solar_new_South > 0) & (solar_new_South <= lambda_S_step_6)) {
				replace lambda_S1 = lambda_S_step_6 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6) & (solar_new_South <= (lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_5 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 ) & (solar_new_South <= (lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_4 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4  ) & (solar_new_South <= (lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_3 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 ) & (solar_new_South <= (lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_2 if TSO == "TenneT_North" 
				}				

				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 + lambda_S_step_2 ) & (solar_new_South <= (lambda_S_step_1 + lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_1 if TSO == "TenneT_North" 
				}
				
				else {
				replace lambda_S1 = lambda_S_step_0 if TSO == "TenneT_North" 
				}
				
				gen z_t = mc_step_2 - lambda_S1 
				sort TSO_num time_utc
				gen lag_z_t = l.z_t			
				
				gen lambda_N_simulated = lambda_S1 + lag_z_t + (`coeff1' - `coeff2') * `deltaK' if TSO == "TenneT_North"	
				gen MB_North_sim = -dAS_dR + lambda_N_simulated + marginal_emissions_EURton if TSO == "TenneT_North"
				
				// Carryforward this variable, so we can use it in calculation for solar_value in South
				gsort time_utc -MB_North_sim
				by time_utc: replace MB_North_sim = MB_North_sim[_n-1] if missing(MB_North_sim)
				
				// Determine possible energy excess in South and adjust value
				gen excess_south = cum_solar_gen - load if  TSO == "TenneT_South"
				gen d_excess_south = excess_south > 0 & !missing(excess_south)
				replace excess_south = 0 if missing(excess_south) 
				
			* CONDITION 1: EXCESS ENERGY + CAPACITY CONSTRAINT NOT BINDING
				gen cond = (excess_south > 0 & excess_south <= `deltaK') 
				gsort time_utc -cond
				by time_utc: replace cond = cond[_n-1] if missing(cond) // within same time interval, carryforward to all TSOs

				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = excess_south * MB_North_sim if TSO == "TenneT_South" & cond == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond == 1	
				
			* CONDITION 2: EXCESS ENERGY + CAPACITY CONSTRAINT BINDING
				gen cond_congest = (excess_south > 0 & excess_south > `deltaK') 
				gsort time_utc -cond_congest
				by time_utc: replace cond_congest = cond_congest[_n-1] if missing(cond_congest) // within same time interval, carryforward to all TSOs
		
				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = `deltaK' * MB_North_sim if TSO == "TenneT_South" & cond_congest == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond_congest == 1	
				
			*** END SECTION ON TRANSMISSION 
			}
	
			else if  ((solar_new > step_6 + step_5 + step_4 + step_3 + step_2) & (solar_new <= (step_1 + step_2 + step_3 + step_4 + step_5 + step_6))) {
				* Apply MB from step 1
				*di "%%%%%%%%%%%%%%%%%%%%%%%%%%"
				*di "In step 1"
				gen MB = -dAS_dR + mc_step_1 + me_step_1
				gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
				gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0 )
				gen MB_MC = cond(cum_solar_gen  < load, mc_step_1 * solar_gen_s, 0)
				gen MB_ME = cond(cum_solar_gen  < load, me_step_1 * solar_gen_s, 0)
				*local in_step_1 `in_step_1' + 1

			*** SECTION ON SIMULATED LAMBDA FOR NORTH - TRANSMISSION **** 
				// NEW: NEED TO DEFINE WHICH lambda_S applies!!
				* While lambda_N = mc6, lambda_S can take on any value from MC0 to MC6: check in loop! * 
				
				** Generate variable to keep track of solar_new_South and carryforward to all TSO 
				gen solar_new_South = solar_new if TSO == "TenneT_South"
				gsort time_utc -solar_new_South
				by time_utc: replace solar_new_South = solar_new_South[_n-1] if missing(solar_new_South) // within same time interval, carryforward to all TSOs
				
				if ((solar_new_South > 0) & (solar_new_South <= lambda_S_step_6)) {
				replace lambda_S1 = lambda_S_step_6 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6) & (solar_new_South <= (lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_5 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 ) & (solar_new_South <= (lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_4 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4  ) & (solar_new_South <= (lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_3 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 ) & (solar_new_South <= (lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_2 if TSO == "TenneT_North" 
				}				

				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 + lambda_S_step_2 ) & (solar_new_South <= (lambda_S_step_1 + lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_1 if TSO == "TenneT_North" 
				}
				
				else {
				replace lambda_S1 = lambda_S_step_0 if TSO == "TenneT_North" 
				}
				
				gen z_t = mc_step_1 - lambda_S1 
				sort TSO_num time_utc
				gen lag_z_t = l.z_t			
				
				gen lambda_N_simulated = lambda_S1 + lag_z_t + (`coeff1' - `coeff2') * `deltaK' if TSO == "TenneT_North"	
				gen MB_North_sim = -dAS_dR + lambda_N_simulated + marginal_emissions_EURton if TSO == "TenneT_North"
				
				// Carryforward this variable, so we can use it in calculation for solar_value in South
				gsort time_utc -MB_North_sim
				by time_utc: replace MB_North_sim = MB_North_sim[_n-1] if missing(MB_North_sim)
				
				// Determine possible energy excess in South and adjust value
				gen excess_south = cum_solar_gen - load if  TSO == "TenneT_South"
				gen d_excess_south = excess_south > 0 & !missing(excess_south)
				replace excess_south = 0 if missing(excess_south) 
				
			* CONDITION 1: EXCESS ENERGY + CAPACITY CONSTRAINT NOT BINDING
				gen cond = (excess_south > 0 & excess_south <= `deltaK') 
				gsort time_utc -cond
				by time_utc: replace cond = cond[_n-1] if missing(cond) // within same time interval, carryforward to all TSOs

				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = excess_south * MB_North_sim if TSO == "TenneT_South" & cond == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond == 1	
				
			* CONDITION 2: EXCESS ENERGY + CAPACITY CONSTRAINT BINDING
				gen cond_congest = (excess_south > 0 & excess_south > `deltaK') 
				gsort time_utc -cond_congest
				by time_utc: replace cond_congest = cond_congest[_n-1] if missing(cond_congest) // within same time interval, carryforward to all TSOs
		
				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = `deltaK' * MB_North_sim if TSO == "TenneT_South" & cond_congest == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond_congest == 1	
				
			*** END SECTION ON TRANSMISSION 	

			}
	
			else {
				* Apply MB from step 0
				*di "%%%%%%%%%%%%%%%%%%%%%%%%%%"
				*di "In step 0"
				gen MB = -dAS_dR + mc_step_0 + me_step_0
				gen MB_tot = cond(cum_solar_gen  < load, MB * solar_gen_s, 0) 			
				gen MB_AS = cond(cum_solar_gen  < load, -dAS_dR * solar_gen_s, 0 )
				gen MB_MC = cond(cum_solar_gen  < load, mc_step_0 * solar_gen_s, 0)
				gen MB_ME = cond(cum_solar_gen  < load, me_step_0 * solar_gen_s, 0)
				*local in_step_0 `in_step_0' + 1
	
			*** SECTION ON SIMULATED LAMBDA FOR NORTH - TRANSMISSION **** 
				// NEW: NEED TO DEFINE WHICH lambda_S applies!!
				* While lambda_N = mc6, lambda_S can take on any value from MC0 to MC6: check in loop! * 
				
				** Generate variable to keep track of solar_new_South and carryforward to all TSO 
				gen solar_new_South = solar_new if TSO == "TenneT_South"
				gsort time_utc -solar_new_South
				by time_utc: replace solar_new_South = solar_new_South[_n-1] if missing(solar_new_South) // within same time interval, carryforward to all TSOs
				
				if ((solar_new_South > 0) & (solar_new_South <= lambda_S_step_6)) {
				replace lambda_S1 = lambda_S_step_6 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6) & (solar_new_South <= (lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_5 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 ) & (solar_new_South <= (lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_4 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4  ) & (solar_new_South <= (lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_3 if TSO == "TenneT_North" 
				}
				
				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 ) & (solar_new_South <= (lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_2 if TSO == "TenneT_North" 
				}				

				else if ((solar_new_South > lambda_S_step_6 + lambda_S_step_5 + lambda_S_step_4 + lambda_S_step_3 + lambda_S_step_2 ) & (solar_new_South <= (lambda_S_step_1 + lambda_S_step_2 + lambda_S_step_3 + lambda_S_step_4 + lambda_S_step_5 + lambda_S_step_6))) {
				replace lambda_S1 = lambda_S_step_1 if TSO == "TenneT_North" 
				}
				
				else {
				replace lambda_S1 = lambda_S_step_0 if TSO == "TenneT_North" 
				}
				
				gen z_t = mc_step_0 - lambda_S1 
				sort TSO_num time_utc
				gen lag_z_t = l.z_t			
				
				gen lambda_N_simulated = lambda_S1 + lag_z_t + (`coeff1' - `coeff2') * `deltaK' if TSO == "TenneT_North"	
				gen MB_North_sim = -dAS_dR + lambda_N_simulated + marginal_emissions_EURton if TSO == "TenneT_North"
				
				// Carryforward this variable, so we can use it in calculation for solar_value in South
				gsort time_utc -MB_North_sim
				by time_utc: replace MB_North_sim = MB_North_sim[_n-1] if missing(MB_North_sim)
				
				// Determine possible energy excess in South and adjust value
				gen excess_south = cum_solar_gen - load if  TSO == "TenneT_South"
				gen d_excess_south = excess_south > 0 & !missing(excess_south)
				replace excess_south = 0 if missing(excess_south) 
				
			* CONDITION 1: EXCESS ENERGY + CAPACITY CONSTRAINT NOT BINDING
				gen cond = (excess_south > 0 & excess_south <= `deltaK')
				gsort time_utc -cond
				by time_utc: replace cond = cond[_n-1] if missing(cond) // within same time interval, carryforward to all TSOs

				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = excess_south * MB_North_sim if TSO == "TenneT_South" & cond == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond == 1	
				
			* CONDITION 2: EXCESS ENERGY + CAPACITY CONSTRAINT BINDING
				gen cond_congest = (excess_south > 0 & excess_south > `deltaK') 
				gsort time_utc -cond_congest
				by time_utc: replace cond_congest = cond_congest[_n-1] if missing(cond_congest) // within same time interval, carryforward to all TSOs
		
				*Adjust values for South (export) and North (apply same MB_sim)
				replace MB_tot = `deltaK' * MB_North_sim if TSO == "TenneT_South" & cond_congest == 1	
				replace MB_tot = cond(cum_solar_gen  < load, MB_North_sim * solar_gen_s, 0) if TSO == "TenneT_North" & cond_congest == 1	
				
			*** END SECTION ON TRANSMISSION 		
			}
	
		} //else
	
	
		** Compare total value across TSOs and update allocation, where installed
		*bysort TSO: egen MB_mean = mean(MB) // total expected MB from current MW of solar
	
		** OBTAIN RANKING IN AUTOMATED FASHION ** 
		egen meanMB = mean(-MB_tot), by(TSO)
		egen MB_order = axis(meanMB TSO), label(TSO)
		qui: tabstat MB_tot , by(MB_order) s(mean sd N) save
		local first   `r(name1)' 
		local second  `r(name2)' 
		local third   `r(name3)' 
		local fourth  `r(name4)'
		local fifth   `r(name5)'
		* su MB if TSO == "`first'"
		di "RANKING TSO: `first', `second', `third', `fourth', `fifth'"
	
	
		// Allocation rule: nest options // 
		** Store total solar installed, cumulative value of solar installed
	
		gen alloc_constraint = 1 if (solar_cap > max_cap_gamma) & TSO == "`first'" 
		replace alloc_constraint = 0 if alloc_constraint!=1 // replaces value in all TSOs
		qui: su alloc_constraint
		local alloc_1 =  `r(mean)'
	
		if `alloc_1' == 0 {
			replace solar_cap = solar_cap - `s' if (TSO == "`second'" | TSO == "`third'" | TSO == "`fourth'" | TSO == "`fifth'")
			* SET back to original value for the TSO's where no allocation - consider constaints!! 
		}
	
		else if `alloc_1' != 0 {
			drop alloc_constraint
			gen alloc_constraint = 1 if  (solar_cap > max_cap_gamma) & TSO == "`second'" 
			replace alloc_constraint = 0 if alloc_constraint!=1
			qui: su alloc_constraint
			local alloc_2 =  `r(mean)'
	
			if `alloc_2' == 0 {
				replace solar_cap = solar_cap - `s' if (TSO == "`first'" | TSO == "`third'" | TSO == "`fourth'"  | TSO == "`fifth'") 
				* SET back to original value for the TSO's where no allocation - consider constaints!! 
			}
	
			else if `alloc_2' != 0 {
				drop alloc_constraint
				gen alloc_constraint = 1 if (solar_cap > max_cap_gamma) & TSO == "`third'" 
				replace alloc_constraint = 0 if alloc_constraint!=1
				qui: su alloc_constraint
				local alloc_3 =  `r(mean)'
	
				if `alloc_3' == 0 {
					replace solar_cap = solar_cap - `s' if (TSO == "`first'" | TSO == "`second'" | TSO == "`fourth'" | TSO == "`fifth'") 
					* SET back to original value for the TSO's where no allocation - consider constraints!! 
				}

				else if `alloc_3' != 0 {
				drop alloc_constraint
				gen alloc_constraint = 1 if (solar_cap > max_cap_gamma) & TSO == "`fourth'" 
				replace alloc_constraint = 0 if alloc_constraint!=1
				qui: su alloc_constraint
				local alloc_4 =  `r(mean)'
	
				if `alloc_4' == 0 {
					replace solar_cap = solar_cap - `s' if (TSO == "`first'" | TSO == "`second'" | TSO == "`third'" | TSO == "`fifth'") 
					* SET back to original value for the TSO's where no allocation - consider constraints!! 
				}	
	
				else {
					// Put in fifth TSO 
					replace solar_cap = solar_cap - `s' if (TSO == "`first'" | TSO == "`second'" | TSO == "`third'" | TSO == "`fourth'" ) 
				} // TSO 5
			} // TSO 4
		} // TSO 3
	} // TSO 2
	
	
		// Store total value and in case "decomposition"
		qui: levelsof TSO, local(T)
		foreach t in `T' {
			replace value_`t' = value_`t' + MB_tot if (solar_cap != solar_capl1) & TSO == "`t'"
			replace value_AS_`t' = value_AS_`t' + MB_AS if (solar_cap != solar_capl1) & TSO == "`t'"
			replace value_MC_`t' = value_MC_`t' + MB_MC if (solar_cap != solar_capl1) & TSO == "`t'"
			replace value_ME_`t' = value_ME_`t' + MB_ME if (solar_cap != solar_capl1) & TSO == "`t'"
		}
		// value_TSO: sums total MB for each block of solar that is reallocated
	
		replace solar_capl1 = solar_cap
		drop dAS_dR alloc_constraint MB_order meanMB MB MB_tot MB_AS MB_MC MB_ME
		drop cond cond_congest excess_south d_excess_south MB_North_sim lambda_N_simulated lag_z_t z_t solar_new_South
	
	} // loop over S (total capacity to reallocate)
	
	
	// AGGREGATE TOTAL BENEFITS 
	qui: su solar_cap if TSO =="50Hertz"
	local solar_cap_50Hertz = r(mean) 
	
	qui: su solar_cap if TSO =="Amprion"
	local solar_cap_Amprion = r(mean)
	
	qui: su solar_cap if TSO =="TenneT_North"
	local solar_cap_TenneT_N = r(mean)

	qui: su solar_cap if TSO =="TenneT_South"
	local solar_cap_TenneT_S = r(mean)
	
	qui: su solar_cap if TSO =="TransnetBW"
	local solar_cap_TransnetBW = r(mean)
	
		
	
	// Value by TSO
	
	preserve
	collapse (sum) value* 
	
	qui: su value_50Hertz
		local value_50Hertz = r(mean)
	qui: su value_Amprion
		local value_Amprion = r(mean)
	qui: su value_TenneT_North
		local value_TenneT_N = r(mean)
	qui: su value_TenneT_South
		local value_TenneT_S = r(mean)
	qui: su  value_TransnetBW
		local value_TransnetBW = r(mean)
		
	// Total of average MB over two year period
	local value_total = `value_50Hertz' + `value_Amprion' + `value_TenneT_N' + `value_TenneT_S' + `value_TransnetBW'	

	qui: su value_AS_50Hertz  
		local MB_AS_50Hertz = r(mean)
	qui: su value_AS_Amprion  
		local MB_AS_Amprion = r(mean)
	qui: su value_AS_TenneT_North
		local MB_AS_TenneT_N = r(mean)
	qui: su value_AS_TenneT_South
		local MB_AS_TenneT_S = r(mean)
	qui: su value_AS_TransnetBW
		local MB_AS_TransnetBW = r(mean)
	
	qui: su value_MC_50Hertz  
		local MB_MC_50Hertz = r(mean)
	qui: su value_MC_Amprion  
		local MB_MC_Amprion = r(mean)
	qui: su value_MC_TenneT_North
		local MB_MC_TenneT_N = r(mean)
	qui: su value_MC_TenneT_South
		local MB_MC_TenneT_S = r(mean)
	qui: su value_MC_TransnetBW
		local MB_MC_TransnetBW = r(mean)		
	
	qui: su value_ME_50Hertz  
		local MB_ME_50Hertz = r(mean)
	qui: su value_ME_Amprion  
		local MB_ME_Amprion = r(mean)
	qui: su value_ME_TenneT_North
		local MB_ME_TenneT_N = r(mean)
	qui: su value_ME_TenneT_South
		local MB_ME_TenneT_S = r(mean)
	qui: su value_ME_TransnetBW
		local MB_ME_TransnetBW	= r(mean)	
	
	egen value_AS = rowtotal(value_AS_*)
	egen value_MC = rowtotal(value_MC_*)
	egen value_ME = rowtotal(value_ME_*)
	
	qui: su value_AS
		local value_AS_tot = r(mean)
	
	qui: su value_MC
		local value_MC_tot = r(mean)
	
	qui: su value_ME
		local value_ME_tot = r(mean)					
	restore	
	
	
	drop lambda_S1
	
	// Fill in Matrice with results
	* Mat 1: Gains
	mat ratios[`number_gammas'*(`m'-1) +`j',1] = `SCC'
	mat ratios[`number_gammas'*(`m'-1) +`j',2] = `gamma'
	mat ratios[`number_gammas'*(`m'-1) +`j',3] = `deltaK'
	mat ratios[`number_gammas'*(`m'-1) +`j',4] = `solarcap_all_50Hertz'
	mat ratios[`number_gammas'*(`m'-1) +`j',5] = `solar_cap_50Hertz' 
	mat ratios[`number_gammas'*(`m'-1) +`j',6] = `solarcap_all_Amprion'
	mat ratios[`number_gammas'*(`m'-1) +`j',7] = `solar_cap_Amprion'
	mat ratios[`number_gammas'*(`m'-1) +`j',8] = `solarcap_all_TenneT_North'
	mat ratios[`number_gammas'*(`m'-1) +`j',9] = `solar_cap_TenneT_N'
	mat ratios[`number_gammas'*(`m'-1) +`j',10] = `solarcap_all_TenneT_South'
	mat ratios[`number_gammas'*(`m'-1) +`j',11] = `solar_cap_TenneT_S'
	mat ratios[`number_gammas'*(`m'-1) +`j',12] = `solarcap_all_TransnetBW'
	mat ratios[`number_gammas'*(`m'-1) +`j',13] = `solar_cap_TransnetBW' 
	mat ratios[`number_gammas'*(`m'-1) +`j',14] = `value_bench_total' / 1e6
	mat ratios[`number_gammas'*(`m'-1) +`j',15] = `value_total' / 1e6
	mat ratios[`number_gammas'*(`m'-1) +`j',16] = `value_50Hertz_bench' / 1e6
	mat ratios[`number_gammas'*(`m'-1) +`j',17] = `value_50Hertz' / 1e6
	mat ratios[`number_gammas'*(`m'-1) +`j',18] = `value_Amprion_bench' / 1e6
	mat ratios[`number_gammas'*(`m'-1) +`j',19] =  `value_Amprion' / 1e6
	mat ratios[`number_gammas'*(`m'-1) +`j',20] =  `value_TenneT_N_bench' / 1e6
	mat ratios[`number_gammas'*(`m'-1) +`j',21] = `value_TenneT_N' / 1e6
	mat ratios[`number_gammas'*(`m'-1) +`j',22] =  `value_TenneT_S_bench' / 1e6
	mat ratios[`number_gammas'*(`m'-1) +`j',23] = `value_TenneT_S' / 1e6	
	mat ratios[`number_gammas'*(`m'-1) +`j',24] = `value_TransnetBW_bench' / 1e6
	mat ratios[`number_gammas'*(`m'-1) +`j',25] = `value_TransnetBW' / 1e6

	
	local ++j
	*di `j'
	
	
	} // LOOP OVER GAMMA

	local ++m
	*di `m'
	
	
} // LOOP OVER deltaK

// Extarct values + plot
mat colnames ratios = SCC gamma deltaK cap_50HZ cap_50HZ_alloc cap_Amprion cap_Amprion_alloc cap_TenneT_N cap_TenneT_N_alloc cap_TenneT_S cap_TenneT_S_alloc cap_Transnet cap_Transnet_alloc v_bench v_alloc /// 
v_50HZ v_50HZ_alloc v_Amprion v_Amprion_alloc v_TenneT_N v_TenneT_N_alloc v_TenneT_S v_TenneT_S_alloc v_Transnet v_Transnet_alloc
svmat ratios, names(col)

		
cap log close 


*
// Extarct main results //
** Calculate additional points in curvature
** For Figure 7: min  0.475, max 0.625, delta 0.05
** Run once for full range and once for additinal data points; adjust file here. 

preserve
drop if missing(SCC)
keep SCC - v_Transnet_alloc
//save "../results/results_interpol.dta", replace
save "../results/results_splitTenneTmain.dta", replace
restore


*******************************************************
// PLOTS 

use  "../results/results_splitTenneTmain.dta", clear // main file
append using "../results/results_interpol.dta" // additional points around .5

sort SCC deltaK gamma 

** Reallocation: SHARE of total capacity
gen cap_tot = cap_50HZ_alloc + cap_Amprion_alloc + cap_TenneT_N_alloc + cap_TenneT_S_alloc + cap_Transnet_alloc // all allocated
gen alpha_50Hertz = cap_50HZ_alloc / cap_tot
gen alpha_Amprion = cap_Amprion_alloc  / cap_tot
gen alpha_TenneT_N = cap_TenneT_N_alloc  / cap_tot
gen alpha_TenneT_S = cap_TenneT_S_alloc  / cap_tot
gen alpha_TransnetBW = cap_Transnet_alloc / cap_tot

label var alpha_50Hertz "50Hertz"
label var alpha_Amprion "Amprion"
label var alpha_TenneT_N "TenneT-N"
label var alpha_TenneT_S "TenneT-S"
label var alpha_TransnetBW "TransnetBW"
label var gamma "{&gamma}"

// Online Appendix Figure D.7 //
foreach K in 0 2000 {
preserve
*keep if gamma <= .7
keep if deltaK ==`K'
tw line alpha_50Hertz gamma, lw(medthick) col(ltblue) || ///
	line alpha_Amprion gamma, lw(medthick) lp(-) col(cyan) || ///
	line alpha_TenneT_N gamma, lw(medthick) lp(.-) col(navy) || ///
	line alpha_TenneT_S gamma, lw(medthick) lp(..--) col(ebblue) || ///
	line alpha_TransnetBW gamma, lw(medthick) lp(_) col(teal) ylabel( , grid) xtitle(`"{fontface "Times":{&gamma}}"') ///
	title("{&Delta} K: `K'") /// 
	ytitle("Share of new solar capacity relative to total solar capacity") ///
	yscale(titlegap(3)) ///
	legend(keygap(*.8) symxsize(*.8))  ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1color) name(alphas_`K', replace)
//graph export "../figures/ratios_alphas_v3_splitTenneT.pdf", replace 
restore
}

graph combine alphas_0 alphas_2000 , rows(1) ///
graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1color) 

graph export "../figures/alphas_v3_splitTenneT.pdf", replace 


** GAINS 
*Total Gains as function of gamma
gen ratios = 100 * ((v_alloc / v_bench) -1)
gen ratios_50Hertz = 100 * ((v_50HZ_alloc /  v_50HZ) - 1)
gen ratios_Amprion = 100 * ((v_Amprion_alloc /  v_Amprion) - 1)
gen ratios_TenneT_N = 100 * ((v_TenneT_N_alloc /  v_TenneT_N) - 1)
gen ratios_TenneT_S = 100 * ((v_TenneT_S_alloc /  v_TenneT_S) - 1)
gen ratios_TransnetBW = 100 * ((v_Transnet_alloc /  v_Transnet) - 1)


// Figure 7 // 
preserve
keep if gamma <= .6
tw line ratios gamma if deltaK == 0 , lw(medthick) col(navy) ||  ///
	line ratios gamma if deltaK == 2000 , lp(-) lw(medthick) col(ebblue) ||  ///
	line ratios gamma if deltaK == 4000 , lp(_) lw(medthick) col(cyan) ||  ///
	line ratios gamma if deltaK == 6000 , lp(..--) lw(medthick) col(ltblue)  ///
	legend(lab(1 "{&Delta} K=0 MW") lab(2 "{&Delta} K=2000 MW") lab(3 "{&Delta} K=4000 MW") lab(4 "{&Delta} K=6000 MW")) ylabel( , grid) ytitle("% of gains""relative to benchmark allocation") xtitle(`"{fontface "Times":{&gamma}}"') ///
	graphregion(lwidth(none) ilwidth(none)) plotregion(lwidth(none) ilwidth(none)) scheme(s1color) name(ratios, replace)
restore
	
	graph export "../figures/ratios_v3_splitTenneT_v1.pdf", replace 
*********************************************************************
